Using the density in QMC

A while ago, I made a small post wondering about Levy constrained search in QMC?
That post was idle speculation, but there have been some papers published seriously exploring the idea for combining DFT and QMC.

The series of papers:
[1] Electronic Energy Functionals: Levy-Lieb principle within the Ground State Path Integral Quantum Monte Carlo (L. Delle Site, L. M. Ghiringhelli, D. Ceperley, May 2012)

[2] Levy-Lieb Principle meets Quantum Monte Carlo (L. Delle Site, Nov 2013)

[3] Levy-Lieb principle: The bridge between the electron density of Density Functional Theory and the wavefunction of Quantum Monte Carlo  (L. Delle Site, Nov 2014)

And two related earlier papers:
Levy-Lieb constrained-search formulation as a minimization of the correlation functional (L. Delle Site, Apr 2007)
The design of Kinetic Functionals for Many-Body Electron Systems : Combining analytical theory with Monte Carlo sampling of electronic configurations (L. M. Ghiringhelli, L. Delle Site, Nov 2007)

(Dates are submission dates to arXiv, not publication dates)

Key issues when combining DFT and QMC:

1. How to compute density functionals from a QMC calculation?
   There are two ways to go about this:
1. Invert the samples to find a functional.
   To understand this we need to know a little bit about Orbital-Free Density Functional Theory.   While the foundation of Density Functional Theory is that the ground state energy is a functional of the density, in practice it is a complicated functional (particularly the kinetic energy piece).   The standard approach is to reintroduce orbitals, which can accurately compute the kinetic energy but also at a cost in performance.   Finding a sufficiently accurate functional without needing orbitals could greatly speed up DFT calculations.  Machine learning techniques are starting to be used to compute complicated functionals (see Finding Density Functionals with Machine Learning )
   This approach could yield useful, transferable functionals.
2. Use the QMC sampling to perform the integral over the functional (see [2]).   There is no need for an explicit expression of the functional.
2. How to use a given density in a QMC calculation?
   In [1], the authors propose the use of Ground State Path Integrals (GSPI), with some additional thoughts on how to sample at fixed density in [2].
   The reason for using GSPI is that the trial wavefunction affects efficiency, but not necessarily the final result, and the addition of a density constraint may speed up sampling.   Compare with a VMC approach, where a minimization of the parameters at the fixed density would still be needed (and of course the accuracy would still be limited by the flexibility of the trial wavefunction).

The final proposal is to make the process iterative - use DFT to compute the density, use GSPI to compute new integrals over the functionals (as in 1B),  use those to recompute the DFT density, etc.

It seems like an intriguing proposal, and I'm curious to see how well it would work on a realistic problem ([2] proposes a helium dimer as a good test case).

Video of SciPy 2011 talk

The video from my SciPy 2011 talk is available here.

At the end of the talk, Andy Terrell asked a question about why this structured approach is necessary - why not just program the steps directly in a script?

I agree that the system I presented does seem very 'heavyweight'. But it feels that doing the transformations directly loses some meta-information that would be useful to keep around, such as which transformations were applied at each step. Additionally, by storing the the pieces of the derivation in data structures, they are potentially easier to manipulate (and emphasizes that they can and should be programmatically manipulated). However, Python's meta-tools are quite powerful, and could probably used to manipulate derivations written directly in a script.


Switching to a different issue, the presented example effectively inlines all the functions into the integral. This is clearly not scalable in the face of larger and more complex expressions. The system needs some way to represent pieces that should transform into a function call in the generated code. I have something working now, but it needs some clean-up work.

Modeling derivations: Questions and Answers

The past two posts (part I and part II) described a system (or style) for scientific programming that I've been working out. This post will have more in a Q&A format (The questions are mostly ones I ask myself - perhaps it should better be called a 'rhetoric & rant' format) An abstract for a talk and paper was accepted to SciPy 2011, so hopefully I can explain it clearly enough by then to get some feedback.

1. What is the target application? Goals?
   
   The target goals are better enabling algorithm research and development (that is, writing a scientific code), and a flexibility to target different runtime systems.
   
   Looking at the examples and case studies on various computer algebra system (CAS) websites, I see many cases where the user modeled some system, solved it, and then used the solution to improve or fix the system. The hard part seems to be modeling the system - this is not my target. I'm interested in problems where the second step is the hard part - and one needs to research new and better solvers.
   

   I'm thinking in terms of a traditional scientific code, and how to automate parts of the development workflow. One big issue is how to use a CAS to assist in the development. Also I'm interested in problems that require high performance computing (HPC) - basically problems that will consume as much computing power and algorithm research as possible. This likely means the code will need to be tuned and altered to support different target systems (multicore workstations, distributed memory parallel clusters, accelerators like GPU's, etc).
   

   
   
2. Why not use a compute algebra system directly? It looks like you're creating a lot of extra work for yourself?
   

   I would like a separation between the specification of an algorithm (for scientific computation, this is largely expressed in equations) and its implementation in the runtime system (along with an automated transformation from the spec to the final system).
   
   Typically a CAS is associated with full featured programming language, but implementing the algorithm in that language doesn't preserve this separation.
   
   As a concrete example, a CAS typically includes facilities for performing numerical integration. What if I want to develop a new numerical integration method - are there any examples on how a CAS can be used to do that?
   
   Another issue with implementing HPC-style problems in a CAS is what happens when the existing solvers in the CAS are too slow, or the CAS system language is too slow (eg, how to scale to large parallel systems). The CAS may do a great deal of work for you, but what happens when you want to move past it? This is the part where code generation comes in.
   

3. What is the role of code generation in this system?
   

   More generally, the 'code generation' step is the 'transformation to the target runtime system' step. The transformation may be trivial, if the final equations can be solved by the CAS facilities, but I'm anticipating the need to ultimately target C/C++ or Fortran.
   

   Also, it seems useful to introduce this method in small pieces into an existing code (converting an entire code would usually be too dramatic of a change to be successful). It would smooth the transition for the generated code to fit the style of the existing code. A template system for generation would be a useful approach. A text-oriented (or basic token oriented system, like the C pre-processor) would work, but they have problems. A system oriented around an actual parse tree of the target language would allow the most powerful transformations.
   
   


4. Is this related to reproducible research in any way?
   

   Packages for reproducible research seem to be focused on the level of capturing parameters that are input to a scientific code, so that the results can be repeated and revisited later.
   Executing a program may be reproducible, but the construction of that program is not.
   
   Aside:
   Reproducibility can be achieved through documentation or automation. The previous standard for scientific work depended primarily on documentation. With computers, automation becomes possible as a working method. (For space reasons, this analysis is very simplistic --- there is more to it.)
   </aside>
   
   This work is intended to make the construction of a code automated. One advantage is in making corrections to mistakes in the derivation - fix it, push a button, and the final code is automatically updated.
   
   This should be a complementary to other systems for reproducible research.

5. Why are you doing this?
   I'm tired of tracking down missing minus signs and misplaced factors of 2 in my code. A significant feature of this approach is that the entire chain of derivation and generation of the program is performed by machine, to eliminate transcription errors. (If this sounds too rosy, rest assured there are plenty of other types of errors to make).
   

Example of modeling derivations

The last post discussed a style of scientific programming that modeled derivations symbolically to create a machine-checked chain from starting equation to final code.

This post will provide an example using the configuration integral for the canonical ensemble in statistical mechanics. I want to evaluate the integral using various methods, starting with grid-based integration methods. Using a grid-based method suffers from the curse of dimensionality, where adding particles raises the time to compute exponentially (which is why Monte Carlo methods are normally used). However, it can still be useful to check Monte Carlo codes with other integration methods for small particle numbers. Also, the free energy (and all associated quantities) is computed directly, where it is more difficult to compute with MC methods. Finally, I'm curious as to how large of a system could be computed on today's hardware.

Sympy is the computer algebra system I'm using. The code for this example is available on github in the derivation_modeling branch, in the prototype subdirectory.

The output (MathML and generated python) can be seen here.

The flow is as follows:

1. Derive the composite trapezoidal rule, starting from the expression for a single interval


2. Generate code for the composite trapezoidal rule


3. Start from the configuration integral (which I've called the partition function in the documents) and transform it by specializing various parameters - 2 particles in 2 spatial dimensions, using the Lennard-Jones interaction potential, etc - until it becomes completely specified numerically.


4. Generate code that evaluates this integral using the generated integration code from step 2

There's still a tremendous amount of work to do on this, but hopefully an outline of this style of programming is visible from the example.

• The MathML pages have a lot of detailed small steps - it would be nice to have a collapsible hierarchy so all the little steps can be hidden.


• Some abstractions should remain during the steps - in particular the potential. The expression for the integrand can get quite messy from explicitly inserting the expression for the potential (and this is a simple case - more complex potentials would make it unreadable).


• More backends for the code generation - C/C++ is the obvious choice in pursuit of higher performance. (The generated python can be run through Shed Skin, for another route to C++). Converting to Theano expressions might be a quick way to try it out on the GPU.

Modeling derivations

A previous post describes how I kept reference versions of important derivations in a document. I would like to move the idea out of LaTeX and into some format more amenable to machine checking. To do this requires some structure for describing and modeling a derivation.

A derivation starts with some base equation, then proceeds through a series of steps transforming it into a form that can be solved numerically (A symbolic/analytic solution would be nice, but not possible in most cases.) The steps can be divided into three categories

1. Exact transformations - rearranging terms, replacing an expression with a different form, the usual steps involved in a proof, etc


2. Approximations - using a truncated series expansion, replacing a derivative with discretized version, etc.


3. Specializations - Specifying the physical or model parameters in order to make the equation numerically solvable - like the number of spatial dimensions, number of particles, interaction potentials, etc.

At the end of the derivation, the result can be evaluated numerically using the CAS (Computer Algebra System - I assume this derivation modeling has been implemented using one) facilities, if simple enough. More complex or time-consuming results (a more likely case) can be further transformed by code generation to C/Fortran (or some other execution system).

One important feature of this workflow is every step in creating the program is checked by machine. This should greatly reduce transcription errors (dropped minus signs and factors of 2 - this is what started the whole project in the first place). Also it moves part of the process of creating scientific programs into the realm of 'ordinary software', and as such makes it subject to normal software engineering techniques - unit testing, etc. Testing scientific (numerical) software is a difficult problem - this should help make it (a little) easier and more reliable.

Some other work I've run across related to describing mathematics:

• Structured derivations were developed for teaching proof in high school mathematics.


• OMDoc is a markup language for describing mathematical structures at a higher level than a single formula (a proof, for instance)

Levy constrained search in VMC?

The Levy constrained search formalism is a useful way to reason about density functional theory. It finds the minimum energy using two nested searches. The inner search seeks the minimum energy wavefunction for a particular density. The outer search seeks the density with the minimum energy.

This could be directly implemented using VMC. The main ingredient is energy minimization subject to a density constraint. Has anyone done this?

It wouldn't be very efficient compared to optimizing the VMC energy directly, but it would make an interesting connection between the methods.

Why scientific computations fail

(I was busy this fall with another project. Now it's time to return to scientific programming)

The previous post covered several levels of tests for assessing the correctness of scientific computation. Now let us look at reasons why the results might not match what is expected.

1. Programming errors
   
2. Mistakes in derivation and transcription of equations
   
3. Violations of assumptions or conventions
   
4. Errors from floating point
   
5. Errors and undesirable behavior from the algorithm (often problem-dependent)
   
   
   
   • Algorithms that have been well-studied by numerical analysts
     
   • Those that have not
     
   
   

Programming errors

The usual sorts of problems in programs (array indices, memory problems, etc). This category of problems can be tested with unit tests and can make use of tools and practices from mainstream software engineering.

Mistakes in derivation and transcription of equations

This is similar to the implementation not matching the specification in traditional software engineering. With scientific computation, though, the "specification" could be more precise in principle (equations and algorithms), and computer algebra systems could aid with performing and verifying steps in the derivation.

Violations of assumptions or conventions

Usually happens during derivation of the equations. A stronger connection between the derivation and the final code might help this.

Errors and instabilities from floating point

Accumulated round-off, loss of precision from subtracting nearly equal numbers, and other problems. For more subtleties than you every wanted to know about, visit W. Kahan's site.
One of the most annoying parts about floating point is lack of good model for understanding (and analyzing) floating point errors. I will return to this issue in a future post.

Errors and instabilities from the algorithm.

This category might be better divided into two subcategories - algorithms that have been well-studied by numerical analysts (like matrix operations, certain ODE's and PDE's, optimization), and those that have not. Many algorithms have been well-studied and have nice packages available to use as black-boxes. However, most scientific programs are a composition of well-known algorithms along with some problem-dependent features as well.


Often times the behavior of the algorithm depends on the details of the physical system (ie, is it a metal or insulator)?


This category also includes science-based and numerical approximations.

Now, how do we go from test case failure to deciding which of these categories the problem falls into (and then, of course, to locating the precise problem)?

• Unit testing will help with the first category. But not so much with the rest. Unit testing is helpful when you want to make sure things haven't changed. The hard part is getting it right the first time.
  
• The different levels of tests in the previous post should give some insight into which category the problem might occur.
  

First I want to look at various models of floating point error and see if we can gain some understanding of the intrinsic uncertainty in a calculation. This seems to be a necessary prerequisite for answering the question of 'how close is close enough?', in order to decide whether a test has actually failed or not.